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Various experiments indicate the existence of a second knee around energy E = 3x 10" eV m 
the cosmic ray spectrum. This feature could be the signature of the end of the galactic component 
and of the emergence of the extragalactic one, provided that the latter cuts off at low energies. 
Recent analytical calculations have shown that this cut-off could be a consequence of the existence 
of extragalactic mag netic fields (Refs. QQl): low energy protons diffuse on extragalactic magnetic 
fields and cannot reach the observer within a given time. We study the infiuence of inhomogeneous 
magnetic fields on the magnetic horizon, using a new semi-analytical propagation code. Our results 
indicate that, at a fixed value of the volume averaged magnetic field {B}, the amplitude of the low 
energy cut-off is mainly controled by the strength of magnetic fields in the voids of the large scale 
structure distribution. 
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I. INTRODUCTION 



Recent developments in our undertanding of the ori- 
gin of cosmic rays have led to the suggestion that the 
cosmic ray spectrum might comprise only two compo- 
nents, one of galactic origin, dominant in the energy 
range E < 10^^ eV, and another of extragalactic origin 
in the energy range E > 10^^ eV In this interpreta- 
tion, the cross-over between the two components marks 
the so-called "second knee" in the all particle cosmic ray 
spectrum (see Ref. Q for a review on cosmic ray data at 
high energies), while the "first knee" at E ^ 2x 10^^ eV 
would be associated with a change of propagation regime 
or the maximal energy of protons at the source. The third 
feature in the cosmic ray spectrum, i.e. the "ankle" at 
E ~ 10^^ eV would be a consequence of pair production 
losses of ultra-high energy protons accelerated in sources 
located at cosmological distances Q. 

On purely phenomenological grounds, this modern 
view is appealing when compared to the more tradi- 
tional interpretation in which the ankle is associated with 
the emergence of an extragalactic cosmic ray compo- 
nent out of a more steeply falling spectrum at energies 
E < 10^^ eV. One of these major advantages certainly 
is the "economy" of sources. It is indeed notoriously 
difficult to accelerate particles beyond E « 10^^ eV in 
supernovae remnants [5|. Therefore, if the first knee cor- 
responds to the maximal energy of protons at the source, 
the fall-off of the galactic component at ~ 10^^ eV would 
naturally be associated with the maximal energy of the 
iron component at the source. On the contrary, if the 
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extragalactic component appears at the ankle, one needs 
to postulate the existence of a third cosmic ray compo- 
nent between ~ 10^^ eV and ~ 10^^ eV (see for instance 
Refs. @ for a recent proposal), or to assume that super- 
novae are able to accelerate particles up to the ankle (see 
for instance Q)- 

The modern interpretation of a transition between the 
galactic and the extragalactic component at the second 
knee does not come without flaws, however. In particu- 
lar, the smooth matching of the galactic and extragalac- 
tic components a.t E ^ 10^^ eV bears the unaesthetic 
look of fine-tuning. Of course, this problem is generic 
to the matching of two distinct components at a point 
where the slope steepens; the introduction of a third cos- 
mic ray component would not help in this respect. As 
far as the galactic component is concerned, the fall-off at 
E > 10^^ eV arises as a direct consequence of the ob- 
servation of the first knee, as mentioned above. However 
one must explain why the extragalactic component van- 
ishes at energies below the second knee. In the original 
scenario of Berezinsky et al. [s!] , this low energy cut-off 
was attributed to physics at the source. In Ref. it was 
suggested to interpret it as the modulation of the extra- 
galactic flux due to a galactic magnetized wind, although 
the calculations of Ref. 9] bring down this cut-off to a 
too low energy, E - lO^''^ _ ^qis eV. 

Another possibility, advocated in Refs. is to relate 
this cut-off with the influence of extragalactic magnetic 
flelds. As demonstrated in these studies, if the inten- 
sity of extragalactic magnetic fields is rather modest, say 
(B) ~ 10"^ G, the diffusion time of particles with en- 
ergy E < 10^^ eV from the closest sources (located at, 
say ~ 50 — 100 Mpc) becomes longer than the age of the 
Universe. This produces a low energy cut-off in the prop- 
agated spectrum at the required location, which allows 
to reproduce a smooth transition at the second knee in 
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agreement with observational data (see [l|). 

It has also been argued that if ultra-high energy cos- 
mic rays comprised a significant fraction of heavy nu- 
clei (> 20%), the scenario of a transition at the second 
knee would loose its merits as the energy losses would 
no longer be able to reproduce the ankle feature [IH ■ 
However even a solar type (or Galactic cosmic ray type) 
chemical composition, with ~ 10 % helium and only 
traces of heavier elements, allows to fit the existing data 
at the ankle with a single powerlaw spectrum at injec- 
tion Q. Furthermore, there is no particular reason to 
expect the source composition to be enriched in metals. 
In any case, future measurements of the chemical com- 
position will tell [Hill. 

Quantizing the influence of extragalactic magnetic 
fields on the spectrum of cosmic rays with energy E ~ 
10^^ eV is not an easy task as the propagation times be- 
come of the order of a Hubble time, hence one must ac- 
count for the effects of expansion. For the sake of sim- 
plicity, Refs. P, i] have thus assumed the magnetic field 
power to be distributed homogeneously in space. How- 
ever, this approximation deserves to be refined since the 
magnetic field is most likely distributed as the charged 
baryonic plasma. Since the scale of inhomogeneity of 
large scale structure in the Universe is comparable to the 
distance to the closest sources, ~ 50 — 100 Mpc, the in- 
homogeneity of the magnetic field may affect the conclu- 
sions of Refs. [JQ- The objective of the present paper 
is precisely to address this issue and to study the sce- 
nario put forward in these references in a more realistic 
extragalactic magnetic field configuration. 

This immediately brings forward the difficulty of defin- 
ing a realistic distribution of large scale extragalactic 
magnetic fields, including the shape and amplitude of 
a turbulent magnetic cascade. From an observational 
point of view, one has been able to measure the strength 
of extragalactic magnetic fields "only" in the core of clus- 
ters of galaxies [l^ ; a bridge of synchrotron emission on 
Mpc scales has been observed in the Coma cluster [l^ . 
Hopefully the SKA project will enlarge considerably the 
dataset on extragalactic magnetic fields [lB| but it is not 
expected to enter operation before 2017. In the mean- 
time, one thus has to rely on theory. Unfortunately, the 
very origin of extragalactic magnetic fields is unknown, 
see Ref. for a review. Furthermore, even if one knew 
exactly the initial conditions that set the configuration 
of magnetic fields at a high redshift, the simulation of 
their evolution throughout cosmic history to the present, 
carrying sufficient accuracy on a large dynamic range of 
spatial scales, remains a formidable task for numerical 
computing. 

In regards of all these uncertainties on the origin of 
extragalactic magnetic fields, on their distribution in the 
present Universe, on the nature and shape of magnetic 
turbulence as well as on the transport properties of par- 
ticles in chaotic magnetic fields, we adopt a simplified 
and parametrized description which allows us to eval- 
uate the effects of the various sources of uncertainties 



on the results. As a by-product of the present study, 
we thus propose a simple and new recipe to build semi- 
realistic magnetic field distributions out of dark matter 
simulations (which can be obtained at a lesser cost than 
MHD numerical simulations) as well as a new transport 
scheme which is more efficient than existing codes in sev- 
eral respects. In particular, it allows to enlarge artifi- 
cally the range of scales on which the magnetic field is 
distributed, hence to model the influence of intergalac- 
tic magnetized turbulence on particle transport. These 
techniques, which are developed in Section [III and in Ap- 
pendix [XI allow us to bracket the possible distributions 
of extragalactic magnetic fields at the present time and 
their impact on the ultra-high energy cosmic ray spec- 
trum. Our results indicate that, at a fixed value of the 
volume averaged (B), the amplitude of the low energy 
cut-off is controlled by the strength Svoid of magnetic 
fields in the voids of the large structure distribution and 
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the source distance scale Us —50 Mpc. The fact that 
our conclusions depend more weakly on other parameters 
characterizing the magnetic field distribution provides an 
adequate a posteriori justification for our semi-analytic 
construction. We also argue that this simulation tech- 
nique offers various advantages over existing full-blown 
MHD simulations of large scale structure formation, at 
least as far as cosmic ray propagation is concerned. 

The paper is organized as follows. In Section [Til we 
present our scheme of inhomogeneous magnetic field sim- 
ulation and the numerical technique of cosmic ray trans- 
port. We compare these techniques to existing simula- 
tions and discuss the advantages and drawbacks of each. 
In Section imi we address the issue of the low energy cut- 
off in various models of extragalactic magnetic fields dis- 
tributions, compute the spectra and compare them to 
experimental data. Section IIVI discusses the limitations 
of our approach and possible future avenues of research. 
Finally, Section FVl summarizes our findings. 



II. PROPAGATION OF HIGH ENERGY 
COSMIC RAYS IN EXTRAGALACTIC 
MAGNETIC FIELDS 

The straightforward way to study the influence of 
extragalactic magnetic flelds boils down to performing 
Monte Carlo simulations of particle propagation in a sim- 
ulated magnetized Universe. This, however, brings in two 
major difficulties, which were alluded to earlier but which 
are rarely discussed in the literature: (i) an accurate nu- 
merical modeling of the transport of charged particles 
in magnetic fields; (ii) an accurate numerical modeling 
of the magnetized volume, including magnetized turbu- 
lence. 

Point (i) deals with the theory of cosmic ray diffu- 
sion, which in spite of a long history and recent major 
progress, has not yet reached a consensus on the trans- 
port of cosmic rays in MHD turbulence (see Ref. [l3| 
for a recent review). Actually, the simulation of parti- 
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cle transport in a well-defined MHD environment is not 
trivial even from a purely numerical point of view. For 
example, Ref. has demonstrated that the interpola- 
tion of the magnetic field from a numerical grid gives an 
erroneous description of particle transport if the Larmor 
radius tl < ^min, where Zmin represents the grid size, 
i.e. the minimum scale of the turbulence inertial range. 
Point (ii) deals with the problem of simulating realistic 
MHD flows on a large range of spatial scales, which also 
constitutes a field of research in its own right. 



A. Magnetic field modeling 

Several pioneering works have studied the propagation 
of cosmic ray s in so-called "realistic" magnetized envi- 
ronments [H, [13, Hlj. These studies have constructed 
the magnetized cube out of a hydrodynamical simula- 
tion of large scale structure formation, which follows the 
magnetic field in a passive way for Ref. (20| . and with 
feedback effects on the matter evolution for Ref. [2ll |. 
The initial conditions for this magnetic field have been 
set at a high redshift (although Ref. [l^l also models the 
production of magnetic fields at accretion shock waves) 
and the overall amplitude of this field has been rescaled 
at the end of the simulation so as to reproduce the ob- 
served strength of magnetic fields in the core of clusters 
of galaxies. This ingenious procedure allows to fix the 
volume averaged magnetic field independently of the ori- 
gin of the magnetic field, although the volume averaged 
magnetic field now depends on the details of the ampli- 
fication process during cluster formation. 

In Ref. [lOl, the authors follow the trajectory of cosmic 
rays using Monte Carlo methods while the authors of 
Ref. ^Tgii derive an upper bound on the typical cosmic ray 
deflection using a semi-analytic transport scheme. Their 
conclusions are radically different: the former authors 
derive a typical deflection of ~ 10 — 20° above 10^" eV 
while the latter find a deflection less than a degree at 
these energies. This discrepancy illustrates the inherent 
complexity of such simulations. The complexity and the 
cost of such numerical simulations are such that it has 
not been possible to elucidate the precise origin of this 
discrepancy yet. It is likely that most of this difference 
is to be attributed to the modeling of the extragalactic 
magnetic field, and to a lesser degree, to the transport 
scheme. 

Figure [H which presents the volume filling factors of 
the magnetic field strength obtained in these numeri- 
cal simulations is particularly instructive (the model of 
Ref. |19| is shown as the dot-dashed line, while the model 
of Ref. '2(5| is given by the long-dashed line) . It reveals 
large differences in the volume averaged magnetic field as 
well as in the spatial distribution of these fields (which 
translates in this figure as a difference in the slopes of 
the volume filling factor) . Again, the origin of this differ- 
ence is not understood. This figure clearly demonstrates 
that the simulated magnetized volumes, despite all the so- 
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FIG. 1: Volume filling factor of the magnetic field in different 
scenarios. In dot-dashed line, the magnetic field simulated by 
Dolag and coauthors in Ref. [l^iiii long-dashed line, that 
simulated by Sigl and coauthors 20]. In solid line, the semi- 
analytic model with B oc p'^^^ , in orange dotted line, B oc 
p. In red dashed line, the model B oc p [l + (p/p)^^]; this 
model simulates a volume with unmagnetized voids. In all 
cases the proportionality factor Bq = 2 nG. 



phistication of the numerical codes used, cannot truly be 
deemed as realistic. It also indicates the need for alterna- 
tive methods to study the transport of high energy cosmic 
rays in extragalactic magnetic fields, in order to provide 
new angles of attack on this difficult problem. This con- 
stitutes one major motivation of the present work, in 
which we develop one such method and apply it to the 
study of the low energy cut-off at energies close to the 
second knee. 

Our magnetized volume is constructed in a simple way 
as compared to Refs. [T^H^li this simplicity offering var- 
ious advantages (and admittedly, several drawbacks) as 
discussed further below. The core of our method is to 
map the magnetic field strength over the gas density us- 
ing an analytical relation B{p) (to be specified later) and 
to distribute randomly the magnetic field orientation in 
cells of coherence length Ic- The gas density itself is 
obtained from a high resolution dark matter numerical 
simulation of large scale structure formation (with stan- 
dard cosmological parameters VIa = 0.7, flm = 0.3 and 
i/o — 70km/s/Mpc). Once the volume has been set up, 
cosmic ray trajectories are simulated as follows. At each 
step, the cosmic ray is supposed to enter a spherical cell 
of coherence of the magnetic field defined by its diameter 
Zc, in which the magnetic field orientation is random. The 
time spent in the cell and the direction of exit of the cos- 
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mic ray are then drawn from semi-analytic distributions 
which simulate the transport of the particle in MHD tur- 
bulence, according to studies carried out in Refs. p^[2^ 
(see Appendix [X] for a detailed discussion) . The particle 
is then moved to another coherent cell and the next step 
is simulated. As we explain in Section IIIIB( we finally 
compute the propagated spectrum in the low energy re- 
gion E < 10^'' eV in a semi-analytic way which allows 
us to model the effect of cosmological expansion over the 
course of propagation from the source to the detector. 
In detail, we use the Monte Carlo simulations of particle 
propagation in the extragalactic magnetic fields at zero 
redshift in order to measure the diffusion coefficients, and 
use existing analytical formulae in order to calculate the 
propagated spectra from the diffusion equation in an ex- 
panding space-time. More details on this latter step are 
provided in Section fill Bl 

Note that the density field of dark matter provides a 
good approximation to that of the gas density in the in- 
tergalactic medium (IGM) on scales larger than a few 
hundred kpc, corresponding to the baryon Jeans length. 
Therefore the overall baryonic gas density field can be 
obtained by smoothing the dark matter distribution by a 
window function of this size. Our dark matter simulation, 
run with the hydrodynamical code RAMSES [11], is 512^, 
with extent 280 Mpc, hence with a grid size ~ 560 kpc: 
in this case the minimum scale of the simulation then 
plays the role of the window function and no smoothing 
is required. Of course, this treatment does not provide 
a perfect description of the gas distribution but we have 
checked that decreasing the resolution by a factor 2 does 
not affect our results. Moreover the gas density field 
serves only as a marker of the magnetic field distribu- 
tion, so that the above error is negligible in comparison 
to the uncertainty surrounding the strength and config- 
uration of the magnetic field. The essential is rather the 
law B{p) which provides the mapping between the mag- 
netized volume and the density field. 

In the case of isotropic collapse, it is well known that 
B cx p^^^ in a plasma of infinite conductivity. This law is 
slightly oversimplistic because it ignores the anisotropy of 
collapse in large scale structures, which results in an en- 
hanced amplification of the magnetic field by shear and 
anisotropic compressive flows |24l |. In detail, the equa- 
tions of ideal MHD lead to the conservation law: 
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which allows to derive the magnification of B from the 
deformation of the density field. If the separation 5q^ be- 
tween two points is mapped into i5q^ — V ■ Sq^ through 
deformation, where V indicates the deformation tensor, 
then B</p< is mapped into I? • B</p<, so that: 
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along one (i.e. collapse on a wall) or two (i.e. collapse 
on a filament) spatial directions. The law B esc p^-^ has 
indeed been observed in the simulations of Dolag and 
co-authors 

Viscosity and shear flows during collapse may also am- 
plify further the magnetic field, leading to departures 
from the law B oc p between regions of very different 
density, in particular the voids and the structures (fila- 
ments and pancakes). For instance both simulations of 
Refs. [H, [131 lead to weaker values of the magnetic fields 
in the voids than would be expected from an extrapola- 
tion of the law B oc p to regions of low density. In order 
to bracket these different effects, we consider several re- 
lations B{p): 
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(3) 
(4) 

(5) 



From Eq. ([^l), it is then easy to derive the law B oc f?^^ 
for isotropic collapse, ox B oc p for anisotropic collapse 



The last model is an ad-hoc modeling of the suppression 
of magnetic fields in the voids of large structure which 
leaves unchanged the distribution in the dense intergalac- 
tic medium (meaning p> (p)). 



B. Why a semi-analytical propagation method? 

At this stage, one should compare the respective mer- 
its and drawbacks of this new method with other existing 
techniques. Concerning the magnetic field distribution, 
our method obviously neglects subtle effects such as the 
amplification of the magnetic field in the vicinity of ac- 
cretion shock waves of large scale structure. However, it 
should be clear that no numerical simulation can claim to 
simulate with accuracy the magnetic field in the vinicity 
of cosmological shock waves due to the intricacy of MHD 
physics at play. The amount of amplification, the coher- 
ence length and the shape of the turbulence spectrum 
remain open questions (see however [25| for a detailed 
discussion of the Weibel instability operating at inter- 
galactic shocks). 

By considering a one dimensional law B{p), our simula- 
tion apparently neglects the influence of the velocity fleld 
on the magnetic fleld ampliflcation. Indeed, B should be 
a multidimensional function that depends on p as well 
as on the velocity field in order to take into account dy- 
namo and shear effects. The three models alluded to 
earlier do actually account for these effects up to some 
extent, as they reproduce the characteristic features ob- 
tained in the numerical simulations that include dynamo 
and shear effects. The choice of a random orientation of 
the magnetic field in each coherence cell also neglects the 
infiuence of large scale motions. Simulations which fol- 
low explicitely the magnetic field indicate that this latter 
tends to be aligned with the principal directions of large 
scale structure, i.e. the axis of the filament for exam- 
ple ^2d\. However, one should note that the simulations of 
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Refs. [19|, |20|, |2g| assume an initial magnetic field with in- 
finite coherence length, so that the final coherence length 
of the final magnetic field along the axis of the filament 
equals the length of the filament (see Fig. 4 of Ref. ^26*| for 
an illustration of this effect). This of course is unrealistic 
(unless one assumes an acausal origin for B) since the 
coherence length of the magnetic field should not exceed 
~ 1 Mpc, the typical turn-around time of an intergalac- 
tic eddy of this size being comparable to the age of the 
Universe [Ulli]. 

In contrast, our simulation presents the advantage of 
simulating this multiple field reversal along the filament. 
The alignment of the magnetic field direction along the 
filament should take place if the coherence length is larger 
than the transverse size of the filament (more exactly the 
typical scale height of the density gradient). If, as is more 
likely, the coherence length is smaller, then the compres- 
sion is similar to planar collapse as far as the magnetic 
field in a cell is concerned, hence the field becomes aligned 
transversely to the density gradient, and not necessarily 
along the filament axis. Our method offers the means 
to include this effect but we leave this investigation to 
future work for simplicity. 

A last point concerning the distribution of the mag- 
netic field in the intergalactic medium is related to its 
origin. The simulations of Refs. [l^ [l^, [2^ set the ini- 
tial conditions at high redshift and ignore other sources 
of magnetic pollution of the intergalactic medium (sec 
Ref. [l6j for a review on the origin of cosmic mag- 
netic fields) , such as galactic outfiows, AGN pollu- 
tion [2^, [30, HH , the amplification of magnetic fields in 
accretion shocks of large scale structure (except Ref. [20| 
which uses a model of Biermann battery effects), tur- 
bulent amplification in the IGM [s^, etc. These sources 
should inffuence the transport of high energy cosmic rays 
in two ways: by modifying the relationship between the 
volume averaged (B) and the value observed in cluster 
cores, and by adding additional scattering centers which 
have been omitted in these simulations. Our simulation 
technique offers the freedom to include such localized pol- 
lution effects in a simple and efficient way: one could in- 
clude these highly magnetized regions in our simulation 
cube by sampling them according to the local matter den- 
sity. For the sake of simplicity, in a first step we ignore 
these additional sources and postpone their study to fur- 
ther work. This choice is conservative in so far as the 
inclusion of localized regions of enhanced magnetic field 
would tend to amplify the magnitude of the low energy 
cut-off, all things being equal. 

Finally, our simulation assumes for simplicity that the 
coherence length is uniform in space, whereas it is 
likely to evolve as a function of the density and veloc- 
ity fields. However, this brings in additional parameters 
which enlarge the parameter space. We believe that at 
this stage, it is more reasonable to study the infiuence of 
Ic by performing different runs with different values of Ic 
and comparing the results. Furthermore, the actual value 
Ic is intimately related to the origin of the magnetic field 



(which sets the initial Ic) as well as to the velocity fields 
which distort the field during the evolution, in particular 
with the upbringing of MHD turbulence. Here as well, it 
should be clear that no simulation can claim to simulate 
these various effects with accuracy. 

The issue of turbulence in the IGM is delicate, be- 
cause the Reynolds number in the intergalactic medium 
may take large or moderate values depending on the en- 
vironment. For turbulent excitation on a length scale L 
at velocity u, this Reynolds number reads |33|: 



Tie ~ 10^ 
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Turbulence is thus probably fully developed in most of 
the IGM, except in the high temperature regions repre- 
sentative of clusters of galaxies. In these regions of high 
kinematic viscosity, the shape and extent of the inertial 
range of turbulence is rather complex, and most likely 
influenced by the strong magnetic field, see Ref. [131 for 
detailed discussions. Turbulence plays a fundamental 
role in the transport of charged particles as well as 
in the reshaping of the distribution of the magnetic 
field, but its incorporation in numerical simulations of 
the gas density is extremely complex. As the largest 
scale of the turbulence cannot exceed a few hundred 
kiloparsecs or a megaparsecs, taking into account just 
one or two decades of inertial range necessitates an 
unrealistically high resolution since the cube size must 
remain larger than the inhomogeneity scale ~ 100 Mpc. 
In this respect, our simulations provide more fiexibility 
because our transport scheme in the simulated magnetic 
field allows to simulate the infiuence of a turbulence 
spectrum down to scales well below tl , see Appendix [X] 

The present simulation technique thus combines sim- 
plicity, efficiency with flexibility, and the approximations 
on which it rests appear reasonable in regards of the un- 
certainties surrounding the origin of extragalactic mag- 
netic flelds, the nature of MHD turbulence and the prop- 
erties of cosmic ray transport in such turbulence. Most 
importantly, its parametrized description allows us to 
test the influence of the various parameters on the re- 
sults in contrast with most other works on this topic. 



III. RESULTS 

The following results were obtained by computing the 
trajectories of 10'^ protons in inhomogeneous magnetic 
flelds mapped according to four models, for many sets of 
energies E, magnetic field characteristic values Bq and 
coherence lengths • We will label in what follows "mod- 
els 1—3" our modeling of B{p) presented in Eqs. ([1H5]). 
We add to these models a last one (model 4) for which 
B cx p^/^ and the level of turbulence i] — {5B'^)/{B'^) <C 
1, where 5B is the inhomogeneous perturbation compo- 
nent of B (defined such as: B = {B) + SB). 
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FIG. 2: Trajectories of protons of different energies (solid 
line: E = lO" eV, dotted line: E = 10^'' eV) in a slice of 
simulated universe. The characteristic magnetic field is taken 
as Bo = 1 nG and the coherence length as Ic — 100 kpc. The 
colorbar on the side indicates the intensity of the magnetic 
field (in log). 



Though Bq and (B) have quite similar numerical val- 
ues, they are not strictly equal (they differ approximately 
by a factor 1.5). (B) represents the volume averaged 
magnetic field and Bq is the proportionality factor in 
models (1—4), so that B = Bq x f{p), where f{p) is 
dimensionless. 

The particles are emitted from 10 different sources cho- 
sen randomly among regions of high baryonic density. A 
detailed description of our code is given in appendix [XI 

We first spot the existence of a magnetic horizon using 
the isotropic collapse magnetic field model (Eq. 2]). We 
then move on to other models, calculate their resulting 
propagation spectra and study the influence of our two 
main parameters: Bq and Ic- 



A. Existence of a magnetic horizon 

Figure [2] shows two examples of proton trajectories 
(sohd hue: E = lO^"^ eV, dotted line: E = 10^^ eV) in a 
slice of simulated universe. The characteristic magnetic 
field is taken as _Bo = 1 nG and the coherence length as 
Ic = 100 kpc. 

Obviously the particles at = 10^^ cV and 
E — 10^^ eV evolve completely differently: the former 
travel in a rectilinear regime without being affected 
by changes in density, whereas the latter experience 
a diffusive propagation. Taking a closer look at the 
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FIG. 3: Root mean square of the distance of 10^ particles 
to their source after one Hubble time {tu ~ 13.9 Gyr) as a 
function of their energy, for Bq — 2 nG and Ic — 300 kpc. 
The solid line represents the root mean square of the distance 
and the surrounding color band its variance. The dotted line 
shows the values obtained from analytical calculations in a ho- 
mogeneous magnetic field (Eq. [9]) and the dot-dashed line the 
threshold energy above which the energy loss time becomes 
< te/2. 



diffusive trajectory, one notices the expected intuitive 
correlation between the fluffiness of the trajectory and 
the clustered regions. 

Figure[3]illustrates these comments in a more quantita- 
tive way. It shows the root mean square of the distance 
of 10^ particles to their source after one Hubble time 
(13.9 Gyr) as a function of their energy, for a character- 
istic magnetic field Bq — 2 nG and a coherence length 
Ic = 300 kpc. Our results are no longer valid beyond 
the dot-dashed line which represents the threshold en- 
ergy above which the energy loss time becomes < tii/2, 
as our simulations do not compute energy losses. 

The first striking remark is that particles of energy 
below E ^ 3 X 10^'^ eV cannot travel farther than a dis- 
tance of a hundred megaparsecs from their sources. This 
corroborates the scenario of Refs. [3] on the existence 
of a magnetic horizon, and extends it to the case of a 
inhomogeneous magnetic field. 

The curve presented in Fig. [3] comprises three distinct 
parts: a diffusive part with a slope of ^ 1/6, a semi- 
diffusive part (slope ~ 1) and a quasi- rectilinear part 
with a slope tending towards zero. These trends can be 
naturally explained by analysing the propagation regimes 
at different energies, in homogeneous magnetic fields. 

The quantity {r^y^'^ plotted in Fig. [3] can be easily 



7 



related to the diffusion coefficient through the equation: 

(r2) = 2Dt^, (7) 

where % is the Hubble time and D the diffusion coeffi- 
cient. This equation follows straightly from the definition 
of D = {Ax"^) /2At, where Aa; represents the displace- 
ment during the time interval At. Thus our computa- 
tion of D in our simulations has a direct influence on the 
shape of the curve observed in Fig. [3l 

As explained in appendix \^ our diffusion coefficient 
is calculated following the results of Casse et al. pjj]. ft 
accounts for both diffusive (rL <C ^c) and semi-diffusive 
(''L > ^c) regimes. In the case of a diffusive regime, equa- 
tion (IA2I) becomes 



(8) 



which corresponds to the standard Kolmogorov diffusion 
regime. Besides, when the Larmor radius tl is somewhat 
greater than the coherence length of the magnetic field, 
we have the well known dependence 



(9) 



Knowing that tl oc EB ^, we get from (|8]) and 



(r2) 



,2\l/2 



OC 



E^/^ at low energies and (r^)^/^ cx E for higher 



energies. 

It is quite surprising that these slopes, expected for 
homogeneous magnetic fields, are also observed in our 
inhomogeneous simulations for the magnetic field model 
and the set of parameters presented in Fig. [31 We will 
show in the following section that this is not true for 
other models and parameters. 

When we get to very high energy {E ~ lO^*'^), the 
slope of (r^)^/^ versus E gets weaker, as particles enter 
the quasi-rectilinear regime. Eq. ([7]) is no longer valid as 
particles never reach the diffusion regime. 



tl and then D using (|A2p . It is then easy to obtain R 
and calculate T using (|B3p . 

For the isotropic collapse model (model 1) and the rep- 
resented parameters {Bq — 2 nG, Ic = 100 kpc), there is 
a noticeable difference between the homogeneous and the 
inhomogeneous cases. The cut-off occurs at lower energy 
for the inhomogeneous case, probably due to voids that 
enable particles to travel farther. 

The previous remark does not stand for a travelled 
distance of 10 Mpc (blue dashed lines). On the contrary, 
the transmission factor in an inhomogeneous magnetic 
field is lower than in the homogeneous case. This is due 
to the influence of the dense environment of the source 
where particles were emitted. On a small scale of 10 Mpc, 
low energy particles have just escaped the high density 
region surrounding their source and cannot propagate as 
far as in the homogeneous case. The influence of the 
environment will be discussed in section HVBI 



B. Calculated spectra 

In order to compare our results with observational 
data, we derive spectra from our simulations in the fol- 
lowing way. Ref. ^ has shown that the solution to the 
diffusion equation in an expanding Universe, assuming a 
constant comoving distance between scattering centers, 
and limiting itself to energy losses by expansion (which 
is correct at energies below ~ 10^*'^ eV) takes the form: 



dE 



(10) 



This solution agrees with Ref. [35|, which derives the 
general diffusion equation in an expanding universe and 



Another illustration of the existence of the magnetic 
horizon is presented in figured) The transmission factor 
is plotted as a function of particle energy for three dis- 
tances to the source (dashed lines: 10 Mpc, solid lines: 
100 Mpc, dotted lines: 300 Mpc). Given an initial source 
position, we propagate protons over one Hubble time. 
At a distance R from the source, we calculate the trans- 
mission factor by taking the ratio between the number 
of particles situated beyond R and the total number of 
particles that were emitted. 

Figure H] clearly indicates the presence of a magnetic 
horizon: for energies below - 2 X lO^"^ eV, only haff of 
the emitted particles reach a distance of 100 Mpc in a 
Hubble time. The cut-off energy is lower than for the 
case represented in Fig. [3] due to the lower value of Ic, as 
will be explained in section UlI Bl 

Thin lines represent the analytical transmission factors 
calculated in appendix [B] using the diffusion coefficient 
implemented in our code (Eq. IA2p , for the homogeneous 
case. For a given energy with a particular set of parame- 
ters, one can calculate the corresponding Larmor radius 
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FIG. 4: Particle transmission factor at various distances from 
the source, as a function of particle energy. Thick lines 
are results from the simulation run with Bo = 2 nG and 
/c = 100 kpc. Thin lines represent the analytical transmis- 
sion factor for a homogeneous magnetic field (Eq. |B3}. 
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c = 100 kpc 



B„ = 0.3 nG 
B„ = 2.0 nG 
B„ = 10.0 nG 



mode 




B„ = 2.0 nG, Ic = 100 kpc 

nnodel 1 : B « p^' 
model 2: B « p"" 
model 3: B « p/(1+p"^) 
model 4: B « f/'^ and r(<<1 



FIG. 5: Upper panel: influence of Bo on the spectra for Ic — 
100 kpc and model 1. Middle panel: influence of Ic on the 
spectra for a fixed value of Bo = 2 nG and for model 1. Lower 
panel: influence of dependence of B over p on spectra (models 
1—4), for fixed values of Bo = 2 nG and Ic — 100 kpc. 
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FIG. 6: Same as Fig. H for Bo = 2 nG and h = 100 kpc, for 
models 1—4. Variances are not represented. 



presents the solutions for various energy losses. 

In Eq. pop above, represents the comoving distance 
to source i, Eg(t, E) the required, energy at time t in 
order to have an energy E at to given the energy losses, 
Q{Eg) the emission rate per source at energy Eg and A 
the comoving "path length" . A is defined as: 



A^ 



dt 



D 



UnEr. 



ait) 



(11) 



where Oo is the scale factor at emission and D the diffu- 
sion coefficient. Physically, A represents the typical dis- 
tance travelled by diffusion, accounting for energy losses. 

In order to calculate A, we first study the dependence 
of D on E, Bq and Ic using our simulations. We find 
that each set of parameters corresponds to a different 
function D{E, Bq, Ic)- We then parametrize the evolution 
of the magnetic configuration as done by Berezinsky & 
Gazizov [36| . as: 

lc{z) = lc{l + z) and Bq{z) ^ Bo{l + z)^-"", (12) 

where m characterizes the MHD amplification of the 
field. For simplicity, we set m to in our calculations. 
This toy model corresponds to a constant comoving dis- 
tance between scattering centers and ignores magnetic 
field amplification during structure formation. In this 
way, we obtain the required dependence of D over t and 
thus the function D[ae Ee/a{t)]. 

The function dEg{t,E)/dE is calculated by inte- 
grating the energy losses, following the calculations 
of Berezinsky et al. [3]. The injection spectrum ex- 
tends from 10^6 eV to £',„ax = 10^° eV. The function 
Q{Eg) — K{Eg/ E^a-x)^^ gives the emission rate per 
source at energy Eg, K being a normalisation factor 
such that / dE EQ{E) — L, with L the total luminosity, 
which is assumed to scale as the cosmic star formation 
rate from Ref. [33|. We will assume in our calculation 
a spectral index of 7 = 2.6. In any case, it should 
be pointed out that the choice of the star formation 
rate has little influence on our spectra, since the effects 
of the magnetic horizon dominates those of the star 
formation history on the low energy part of the spectrum. 

At higher energies, when the comoving light cone dis- 
tance r{t) = J^°dt'/a{t') becomes smaller than \{t,E), 
the propagation is no longer diffusive and enters the rec- 
tilinear regime. In this case, the propagated spectrum is 
given by: 



Jicct{E) — - — '^'^ 



1 1 dEgit„E) 



An ^ Anrf 1 



dE 



Q{Eg{t,,E)), 



(13) 



where ti is related to by — J*° dt' /a{t'), ri and Zi 
denoting the comoving distance and redshift of the i^^ 
source. The factor 1/(1 -I- Zi) was omitted in Ref. [l|, 
but it has no influence whatsoever as z <C 1 when the 
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Energy [eV] 

FIG. 7: Total spectra (galactic + extragalactic) compared to data. Each row corresponds to a model and a set of parameters. 
Caution: Hs = 10"" Mpc"^ for the third row. The left panels show KASKADE, Akeno and AGASA data. The right panels 
show KASKADE, HiRes-1, HiRes-2 and Fly's Eyes data. Solid lines represents the median values of the total flux, dot-dashed 
lines the separate galactic and extragalactic components and the dotted lines the upper 75"^ and lower 25**^ percentiles for the 
magnetic cut-off of the extragalactic flux. 



rectilinear regime is reached. 

Figure [H] presents the influence of the i3(p)-models and 
of parameters Bq and Ic on the magnetic cut-off. Only 
the diffusive part of the spectra is represented there and 
the fall-off of the curve around ~ 10^^ eV corresponds 
to the transition between the diffusive and rectilinear 
propagation regimes. We assume continuously emitting 
sources with density Us = 10~^ Mpc~'^ and plot the 
median spectrum obtained over 100 realisations of the 
source locations. For each realisation, the location of 



the first hundred sources were uniformly sampled. For 
farther sources, the continuous source approximation is 
valid and it was used numerically. 

The upper panel shows the intuitive result that the 
greater the mean magnetic field, the steeper the cut-off. 
Of course this law is not restricted to model 1 but is also 
valid for models 2—4. 

The middle panel shows interesting features that are 
in agreement with Eqs. ([5]) and For a fixed value of 
-Bo, for low energies, particles are in the diffusive regime 
[see Eq. ([5])] and (r^) scales with the coherence length 
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as Zc ■ For higher energies particles are in the semi- 
difFusive regime [see Eq. and (r^) (x In other 

words, the spectrum cuts off more steeply for lower values 
of Ic for low energies and for greater values of Ic for high 
energies. This can be seen on figure [Sj for 1^ — 300 kpc, 
the spectrum cuts off at high energy but the slope is 
shallow for low energies, whereas for Ic = 30 kpc, the 
slope is steep at low energies but the cut-off starts at 
lower energies. 

In view of these trends, one will have to find a good 
compromise in order to obtain satisfactory fits to the ob- 
servational spectra. 

The lower panel illustrates the shape of the cut-off 
for the four models previously described. Models 3 and 
4 present a much shallower slope compared to models 
1 and 2. The almost total absence of magnetic field 
in the large scale structure voids for model 3 and the 
cancellation of turbulence in model 4 can explain this. 
We also notice that model 1 which has a higher magnetic 
field intensity in the voids cuts off in a steeper way than 
for model 2. 

For a better understanding of the trends seen in the 
lower panel of Fig. [5l we plot in figure [6] the root mean 
square of the distance of 10'^ particles to their source 
after one Hubble time as a function of their energy, as in 
Fig. [31 for models 1—4. Variances have the same width 
for all models; we did not represent them for clarity. 
As already mentioned, the functions represented in this 
figure are closely related to the diffusion coefficient D 
[see Eq. ([7])], which is required to calculate the spectra 
in the diffusive regime. 

In figure [7] we present the total spectra (galac- 
tic-|-extragalactic) compared to the data, for our param- 
eter fit for each model. As for figure [5l we draw the 
median spectrum (dot-dashed line) obtained over 100 re- 
alisations of the source locations. The upper and lower 
dotted curves show the 75'^ and 25"^ percentiles around 
this prediction, meaning that only 25% of spectra are 
higher or lower respectively than indicated by the curves. 
This uncertainty is related to the location of the closest 
sources. As explained in Rcf. [l|, we draw a straight 
dashed line in the region slightly above 10^* eV, where 
the propagation is neither rectilinear nor diffusive (see 
Ref. [H for more details on this transition zone). 

The galactic cosmic ray component is modeled as fol- 
lows. Supernovae are accepted as standard acceleration 
sites, yet it is notoriously difficult to explain acceler- 
ation up to maximal energy 10^^ eV. Thus it is as- 
sumed that the knee sets the maximal acceleration en- 
ergy for galactic cosmic rays: in this conservative model, 
the spectrum of species i with charge Z takes the form 
jziE) ~ {E/Ez)-''- exp{-E/Ez), with 7, - 2.4-2.7, a 
species dependent spectral index, Ez ~ Z x 2 • 10^^ eV 
[38| . The total galactic component is obtained as the sum 
of elemental spectra, each adjusted to KASKADE data 
as described in Ref [IJ. 



We use the data of six major experiments that 
measured the cosmic ray fluxes in our regions of interest: 
KASKADE (2004 data), with an energy range going 
from 10^"^ to 10^^ eV [38|, Akeno from lO^^ to lO^^ *^ eV 
,3a], AGASA from 10^^-^ to lO^^-^ eV 40], HiRes I and II 
from 101^-3 to 10^0 eV [U and Fly's Eyes from IQi^-^ to 
10^° eV [12]. We split these data in two sets in order to 
account for the discrepancy between HiRes and AGASA. 
This enables us to have two different normalizations for 
the extragalactic flux on the left and right panels. The 
normalization of KASKADE data remains the same for 
both sets. 



Four main points emerge from Fig. [71 (i) The sec- 
ond knee feature appears more or less clearly in the four 
models, but ultimately remains quite robust to model 
changes, (ii) However, again, the influence of the mag- 
netic field intensity in voids is obvious: even with a source 
density of rig = 10~® Mpc^"^, the goodness of fit of model 
3 with the observed spectra is only marginal. This situa- 
tion is clearly improved in the other models, especially if 
we consider the uncertainty on the position of the closest 
sources. (Hi) One might also notice that this last element 
has a considerable impact on the cut-off energy, much 
more than in the case of the homogeneous magnetic field 
of Ref. [l|. This is due to the presence of the diffusive 
regime at the low energy tail. One can indeed observe 
in Fig. [51 the flat diffusive locus at low energies for mod- 
els 1 and 4. Phenomenologically, one understands that 
for these models, a slight change in the closest source 
distance can influence greatly the flux of low energy par- 
ticles, (iv) Finally, comparing our plots for AGASA and 
HiRes data, we conclude that the fits are better for the 
latter. The higher slope above the second knee break 
point in the HiRes data as well as the gap of data be- 
tween the KASCADE and HiRes ranges make the fitting 
easier. 



One should emphasize, however, that the above fits 
were obtained by hand, not by any optimization proce- 
dure due to the computing time required to compute one 
spectrum. Therefore, the spectra shown above do not 
strictly speaking represent the best fit to the data. Fur- 
thermore, one should also exert some caution when com- 
paring datasets from different experiments. In Fig. [71 
we chose to plot separately the AGASA and HiRes data 
because of the well-known discrepancy, but one cannot 
exclude a discrepancy between the energy scales of KAS- 
CADE and HiRes for instance, which would shift one 
dataset with respect to the other. Given all these un- 
certainties, the fits shown in Fig. [71 appear satisfactory, 
except maybe that of model 3 in which the cut-off always 
appears too mild. 
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IV. DISCUSSION 
A. Current limitations 

As already discussed in section [TTl our simulations do 
not take into account several different features of extra- 
galactic magnetic fields, both for the sake of simplic- 
ity and because they are in any case poorly understood 
and poorly constrained. We thus mentioned that our 
fields are related to the gas distribution according to 
three models (Eqs.HHl]), which are one-dimensional, that 
no magnetic source is included, and that the coherence 
length is assumed to be uniform in space. 

Another point that should be underlined is that we 
propagate our particles in a static universe, represented 
by the final output (at z = 0) of a cosmological simula- 
tion. In other words, the magnetic fields do not evolve 
in time during our simulations. The universe being more 
dilute at higher redshifts, the effects of inhomogeneous 
magnetic fields may be less important. A way of im- 
proving our results could be to propagate directly parti- 
cles in an evolving magnetic field. Such a method would 
however be very time consuming and as explained in sec- 
tionini subject to too many uncertainties. One could also 
apply our semi-analytical propagation method to series 
of snapshots of the density of the universe at various z. 
But again, one stumbles over our lack of knowledge about 
extragalactic magnetic fields: we have no hint of the evo- 
lution in time of the relation B(^p). Considering all these 
uncertainties, our restriction to a simple static universe 
thus appears reasonable. 

Finally, we did not account for energy losses during the 
simulations which serve us to "measure" the diffusion 
coefficient, but included them in the calculation of the 
spectra presented in Section IIII Bl For energies below 
E ^ 10^* eV that are of interest to us in this paper, 
only expansion losses play a noticeable role (see (3|), at 
least at low redshift [z < 1). Above ^ 10^^ eV and 
for greater redshifts, energy losses by photopion and pair 
production are no longer negligible. Accounting for these 
energy losses should soften the low energy part of Fig. [3] if 
the magnetic field does not evolve strongly with redshift. 
Indeed, some of the particles at E ^ 10^^"^^ eV actually 
result from higher energy particles that lost their energy. 
The greater distance travelled by these particles before 
losing their energy would tend to raise the rms of the 
distance to the source for low energy particles. Let us 
stress again that in order to model this effect, one would 
need to follow as well the evolution of the magnetic field 
with redshift. 



B. Signatures 

Firstly, we showed in section IIII Bl that only some 
particular types of magnetic fields were able to reproduce 
the data, in the context of our study. Namely, for a 
source density of rig = 10~^ Mpc~'^, the voids of large 
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FIG. 8: Mean deflection angles at various distances from the 
source, as a function of particle energy, for model 1. Solid lines 
are the results of our simulation with Bq = 2 nG, Ic ~ 300 kpc 
and using model 1 (Eq. |4|). Dashed lines are the analytical 
values calculated in Ref. [43| (Eq. [14)1 . 



scale structures should have a certain level of magneti- 
sation, and (B) should be roughly comprised between 
0.3 and 10 nG, 2 nG being a overall satisfactory average 
intensity. These numbers should be taken cautiously, 
remembering all the limitations and unknowns that 
affect these kind of simulations, as stated in the first two 
sections. 

We calculated the Faraday rotation measure (RM) for 
our four magnetic field models with a characteristic mag- 
netic field of Sq = 2 nG. Having sampled 10^ lines of sight 
in our simulation cube, we calculated the median of the 
RMs along them. For our models, the power laws of 
median(RM) versus the distance are steeper (slope ~ 1) 
than that expected for a homogeneous magnetic field, 
for which the integration of RM is equivalent to a simple 
random walk (slope ~ 1/2). 

We find that at a cosmological distance of 1 Gpc, the 
median of our RMs is of order ~ 0.03 rad/m^ for model 1 
and of ^ 0.1 rad/m^ for models 2 and 3. These values 
are consistent with the current observations of RMs that 
predict an upper limit of 5 rad/m^ p^ . It should be 
remarked however that the RMs calculated here are sub- 
ject to high variations according to the concentration of 
matter along the line of sight. Though the distribution 
of the RMs is sharply peaked around 0, with most of the 
RMs in the narrow interval of [—0.5, 0.5] rad/m^, we still 
find some punctual cases where the RM can diverge from 
20 up to 2000 rad/m2. 

The use of a median value of RM enables us to get rid 
of the undesirable lines of sight that cross high density 
clusters and that induce these divergences. These few 
lines of sight have a dominant contribution especially if 
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we calculate the variance or the root mean square of RM, 
leading to very high artificial values. 

Note that our rotation measures are again calculated 
for magnetic fields that do not evolve in time. Hence 
our median values can be considered as upper limits, as 
far as relatively low density regions are observed. 

Mean particle deflection angles induced by the mag- 
netic fields of model 1 are presented in Fig. [5] At a 
given distance from the source, we calculate the deflec- 
tion angle between the arrival direction and the line of 
sight to the source. We stop computing the angles when 
the energy loss distance becomes greater than the linear 
travelled distance. 

Our curves compare quite well to the analytical de- 
flection angles calculated by Waxman and Miralda- 
Escude (dashed hues): 




where r is the distance to the source. For all distances, 
the curves deviate from the analytical model at low en- 
ergies, when diffusion becomes important, and saturate 
at 90°. 

The deflections obtained for cosmological distances 
at high energy are quite moderate for model 1. We 
calculated that it is also the case for models 2 and 4 
(deflections are slightly amplified in model 3). For a 
particle energy of 5 x 10^^ eV and a magnetic field of 
Bq = 2 nG and Ic — 300 kpc, we find that the defiection 
is of order ^ 3 — 5° at 100 Mpc for models 1 and 2, 
and of ~ 8° for model 3. These results are consistent 
with the observations of doublets and triplets of events 
by recent experiments and leave room for doing cosmic 
ray astronomy. Our models would thus be in agreement 
with the detection of counterparts at energies around 
the GZK cut-ofl[. 

Recently, a study related to the present work appeared, 
claiming that partial confinement in magnetic fields sur- 
rounding the source plays an important role in the cut-off 
at low energy [44| . This possibility had been put forward 
in Ref. [1], where it was further shown that the time 
of escape from the dense source environment could be 
non negligible only in a rather contrived situation, since 
it requires B > lfiG{lc/10 kpc)!/^ (L/lOO kpc)-\ with 
L the characteristic scale of the magnetic field spread 
around the source. 

We do not find such a strong effect in our present sim- 
ulations, where the time of confinement remains of order 
~ 1 Gyr for a particle of energy 10^^ eV, or ~ 300 Myr 
for a particle of energy 10^^ eV. This effect obviously de- 
pends largely on the source environment, and on the loca- 
tion of the source. Ref. [3] samples the source locations 
according to the baryon density and therefore tends to 



favor high density (cluster) regions. Since the magnetic 
field in the simulations of Ref. [4J| is already quite strong 
(see Fig.[T]), this explains the magnitude of the effect. As 
noted in Ref. [T], the search for counterparts will allow 
to confirm or exclude this effect by studying the environ- 
ment of the sources. 

Ref. jl3| also calculates spectra for inhomogeneous 
magnetic fields. However, unlike in our work, Ref. (4^ 
does not account for the evolution of the magnetic field 
due to expansion during the propagation. 

Nevertheless it is interesting that Ref. [44] finds the 
spectrum to maintain its 'universal' shape in the region 
of moderate energies, where the transition between the 
diffusive and the rectilinear regimes occurs. This is one 
region which we cannot probe using the semi-analytical 
technique of spectrum reconstruction used in the present 
work; the results of Ref. [44] justify our interpolation of 
the spectrum in this region. 



V. CONCLUSION 

We developed a new method combining an efficient 
propagation scheme and a simple recipe to build semi- 
realistic magnetic field distributions. We map the mag- 
netic field following the baryon density distribution ac- 
cording to three different models and propagate particles 
from cell to cell, taking into account the inner turbulence 
of each cell, as well as its global magnetic intensity. This 
method is much faster than classical trajectory integra- 
tions. 

Under the assumption that the emergence of the extra- 
galactic component occurs at the second knee, we demon- 
strated that it was possible give rough limits for some 
key parameters {{B), Ic), by studying their effects on the 
magnetic horizon. 

For our models assuming isotropic or anisotropic col- 
lapse, with or without turbulence (models 1, 2 and 4 
described in sections [H] and , we find that our cal- 
culated spectra fit the data satisfactorily. Numerically, 
for a source density of ris = 10~^ cm"'' we find that an 
average magnetic field (B) = 2 nG is a reasonable value 
for the three models cited above, and coherence lengths 
of 100 kpc (for models 2 and 4) up to 300 kpc (model 1) 
provide a good agreement with the data. These num- 
bers should still be taken cautiously, remembering the 
limitations discussed throughout this paper. 

We showed that the validity of this scenario depends 
on other parameters (relative normalisation of data sets, 
source density) but eventually, the strongest constraint 
comes from the rate of magnetic enrichment of the low 
density intergalactic medium (voids). We saw indeed 
that model 3, which simulates a volume with unmag- 
netized voids has a marginal goodness of fit with the 
observed spectra, even with a low source density. 

Ultimately, therefore, the success of this scenario for 
the transition between the Galactic and extragalactic 
cosmic ray components depends on the very origin of 
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intergalactic magnetic fields, and on whether the voids 
of large-scale structures have remained pristine or not. 
Interestingly, this question is related to the ongoing de- 
bate on the enrichment of the underdense intergalactic 
medium in metals, since galactic winds carrying metals 
also carry significant magnetic fields. Detailed studies of 
the intergalactic medium as well as progress on extra- 
galactic magnetic fields in the coming decade will shed 
light on this issue. 
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APPENDIX A: NUMERICAL TECHNIQUES 
1. Setting up the magnetic field 

Our basic assumption is that the magnetic field and 
the matter density have a similar spatial distribution at 
large scales. For the reasons detailed in section [TTl we 
map the magnetic field following four models (Eqs.HHS]). 
The parameter Bq mentioned along this paper refers to 
the proportionality factor of Eqs. (|3H5]) and indicates the 
caracteristic intensity of the field. 

We calculate our magnetic field by applying these for- 
mulae to a three dimensional dark matter overdensity 
map (at redshift z = 0) generated by the hydrodynam- 
ical code RAMSES [23|. This cosmological simulation 
was based on a the ACDM model that assumes a flat, 
low density universe, with fim = 0.3, I^a = 0.7 and Hub- 
ble constant h = i/o/(100 km/s/Mpc)= 0.7. It models a 
200 Mpc./i^^ comoving periodic cube split in 512'^ cells, 
where the dark matter overdensity is computed. We do 
not resolve structures below Jeans length, which implies 
that we can identify the computed dark matter distribu- 
tion to a gas distribution. 

While propagating our particles, the magnetic field at 
a given position is computed using the overdensity of the 
nearest grid point. 

2. Particle propagation 



of our particles. The initial direction of the impulsion of 
a particle is also drawn randomly. 

a. Computing the trajectory 

Solving the equation of motion for each particle dur- 
ing an entire Hubble time using a Runge-Kutta method 
can be very time consuming. An alternative method con- 
sists in assuming that a proton travels through adjacent 
spheres of diameter Ic, in which the magnetic field has 
a certain level of turbulence (see Fig. 0. We calculate 
analytically the time of escape from each sphere and sam- 
ple the deflection angle of the particle after each sphere 
from a normal law where the mean deflection m and its 
variance s depend on the Larmor radius of the particle. 
This sphere-crossing method is much faster than a direct 
integration of the trajectory. It also enables us to take 
into account the low level turbulence for scales smaller 
than ^c- 

A particle can cross a sphere in two extreme ways: 
either p = ^l/^c 3> 1 and the particle goes nearly straight 
through the sphere, or p <^ 1 and it wanders in the sphere 
for some time. In this latter case, the time spent in the 
sphere ti ~ Ic/c. 

In the Kolmogorov regime in which the diffusion length 

1/3 2/3 

^scatt ~ Ic with tl < /c, one always have /scatt < 
Ic hence the particle enters the diffusive regime in the 
sphere before exiting. The time spent in the sphere then 
depends on the diffusion coefficient D. 

On average, a particle diffusing through the sphere 
travels a linear distance Ic/V^, so that the time of es- 
cape reads: 



D is computed according to the results of Casse et 
al. ,18-5, who performed Monte Carlo simulations of par- 
ticle propagation in stochastic magnetic fields to measure 
the spatial diffusion coefficients. Their data for full tur- 
bulence can be fit by the approximate relation: 

We use this formula in our code to calculate t. 



For each set of parameters, we propagate 10^ protons 
emitted from 10 different sources. The positions of the 
sources are chosen as follows: we first select the grid 
points in our cube that have an overdensity p > 10. 
Smoothed over our grid of 512'^ cells for 200"^ Mpc'^, over- 
densities of ~ 10^ (minimum density of galaxies) would 
indeed correspond to overdensities of ~ 10. Thus the 
selected regions have a good probability of belonging 
to massive halos. We randomly choose ten grid points 
among this subset and label them as the initial positions 



6. Deflection angle calculation 

Our spheres are located so that the particle always 
enter radially. The defiection angle due to the crossing 
of a sphere is calculated with respect to this entering 
direction (see Fig. [S]). 

The cosine of the deflection angle cos 9 varies accord- 
ing to the rigidity p = tl/^c- For p <ti 1 the particle 
is in a diffusive regime, which implies that cos9 has a 
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100 kpc Bo = 2 nG 



FIG. 9: Sketch of particle propagation in magnetic fields as 
modeled in our simulations. Particles enter spheres of di- 
ameter Ic along a radial direction and escape at calculated 
positions (marked with red dots). A global magnetic field in- 
tensity, that can be related to the largest scale of turbulence, 
is associated to each sphere (labeled Bi, with ,63 > B4,). Tur- 
bulence is taken in account inside each sphere through the 
calculation of the escape time and position. In this figure, 
cos 83 is sampled from a normal law with a large variance and 
cos with a smaller variance. 



uniform distribution over [—1, 1]. For the other extreme 
case, namely p 3> 1, cos^ can be sampled from a normal 
law of mean m and variance s. 

In order to identify the functions m = m{p) and s — 
s{p) in the quasi-rectilinear regime, we integrated the 
trajectories of 10'* particles in a sphere, using a Runge- 
Kutta method, for values of p ranging from 2 to 100. This 
gives us a trend for the mean deflection angle {6) and its 
variance S9. 

Then, in our simulations, we use the calculated prob- 
ability law given by m{p) and s{p) in order to draw the 
direction of exit, and we move the clock forward by the 
time spent in the sphere. 



c. Checking the validity of our code in a homogeneous case 

Our code was checked on a broad range of parame- 
ters in the analytically calculable homogeneous case, i.e. 
we compared the numerical simulations using the above 
sphere crossing trajectories and a magnetic field strength 
equal in each sphere with the analytical results for full 
turbulence and homogeneous magnetic power. We did 
obtain through the numerical calculations the expected 
slopes of 1/6 and 1 for Fig. [3] as predicted by the analyt- 
ical calculations. The transmission factors also fit nearly 
exactly those calculated in Eq. (jB3|) . see Fig. [TOl 




1.0 

R/(4Dt)'/' 



10.0 



FIG. 10: Fraction of particles located at a linear distance 
greater than R after one Hubble time in a homogeneous mag- 
netic field, as a function of i? = R/y/ADt as defined in (|B4|l . 
for various rigidities p — ri^jlc- Solid lines are results from 
our simulations and dotted lines the analytical expression cal- 
culated in appendix IB] 



given time i, in a diffusive propagation regime and in a 
homogeneous space. The density of particles at a position 
r and at time t can be written: 



n(r, t) 



(47rL>i)3/2 



exp 



4L>i 



(Bl) 



where D is the diffusion coefficient. The fraction of parti- 
cles beyond R a,t t can thus be obtained by the following 
integral: 



47rr^ dr ■ 



1 



exp 



Integrating by parts, we obtain: 



T=y^i?e-* + erfc^ 



ADt 



where 



erfc(a::) — —= 



and R = 



R 



(B2) 



(B3) 



(B4) 



APPENDIX B: TRANSMISSION FACTOR FOR 
DIFFUSION IN HOMOGENEOUS SPACE 



We calculate in this section the fraction of particles 
that are located at a linear distance greater than i? at a 



15 



[1] M. Lcmoinc, Phys. Rev. D 71, 083007 (2005). 

[2] R. Aloisio and V. Berezinsky, ApJ 625, 249 (2005). 

[3] V. Berezinsky, A. Gazizov, and S. Grigorieva, Phys. Rev. 

D 74, 043005 (2006), arXiv:hep-ph/0204357. 
[4] M. Nagano and A. Watson, Rev. Mod. Phys. 72, 689 

(2000). 

[5] P.-O. Lagago and C. Ccsarsky, A & A 125, L249 (1983). 
[6] R. Budnik, B. Katz, A. MacFadyen, and E. Waxman 

(2007), arXiv:astro-ph/0705.0041. 
[7] A. R. Bell and S. G. Lucek, MNRAS 321, 433 (2001). 
[8] H. J. Volk and V. N. Zikashvilii, A & A 417, 807 (2004). 
[9] H. Muraishi, S. Yanagita, and T. Yoshida, Prog. Th. 
Phys. 113, 721 (2005). 
[10] V. S. Berezinsky, S. I. Grigorieva, and B. I. Hnatyk, As- 

tropart. Phys. 21, 617 (2004). 
[11] D. AUard, E. Parizot, A. V. Olinto, E. Khan, and 

S. Goriely, A & A 443, L29 (2005). 
[12] D. AUard, E. Parizot, and A. V. Olinto (2007), 

arXiv:astro-ph/0703633. 
[13] P. P. Kronberg, Rep. Prog. Phys. 57, 325 (1994). 
[14] K. T. Kim, P. P. Kronberg, G. Giovannini, and T. Ven- 

turi. Nature 341, 720 (1989). 
[15] B. M. Gaensler, R. Beck, and L. Feretti, New Astr. Rev. 

48, 1003 (2004). 
[16] L. M. Widrow, Rev. Mod. Phys. 74, 775 (2002). 
[17] A. Lazarian, J. Cho, and H. Yan (2002), arXiv:astro- 
ph/0211031. 

[18] F. Casse, M. Lemoine, and G. Pelletier, Phys. Rev. D 65, 
023002 (2002). 

[19] K. Dolag, D. Grasso, V. Springel, and I. Tkachev, JKAS 
37, 427 (2004). 

[20] G. Sigl, F. Miniati, and T. A. Enfilin, Phys. Rev. D 70, 
043007 (2004). 

[21] K. Dolag, D. Grasso, V. Springel, and I. Tkachev, Jour- 
nal of Cosmology and Astro-Particle Physics 1, 9 (2005), 
arXiv:astro-ph/0410419. 

[22] J. Candia and E. Roulet, JCAP 0410, 007 (2004). 

[23] R. Teyssier, A&A 385, 337 (2002). 



[24] M. Bruni, R. Maarteens, and C. G. Tsagas, MNRAS 338, 
785 (2003). 

[25] R. SchUckeiser and P. K. Shukla, ApJ 599, L57 (2003). 
[26] M. Briiggen, M. Ruszkowski, A. Simionescu, M. Hoeft, 

and C. Dalla Vecchia, ApJ 631, L21 (2005). 
[27] E. Waxman and J. Bahcall, Phys. Rev. D 59, 023002 

(1999). 

[28] R. Aloisio and V. Berezinsky, ApJ 612, 900 (2004). 

[29] M. J. Rees and G. Setti, Nature 219, 127 (1968). 

[30] S. R. Furlanetto and A. Loeb, ApJ 556, 619 (2001). 

[31] Gopal-Krishna and P. J. Wiita, ApJ 560, L115 (2001), 
arXiv:astro-ph/0108117. 

[32] R. M. Kulsrud, R. Cen, J. P. Ostrikcr, and D. Ryu, ApJ 
480, 481 (1997), arXiv:astro-ph/9607141. 

[33] K. R. Lang, Astrophysical formulae (Astrophysical for- 
mulae / K.R. Lang. New York : Springer, 1999. (Astron- 
omy and astrophysics library, ISSN0941-7834), 1999). 

[34] A. A. Schekochihin and S. C. Cowley, Physics of Plasmas 
13, 6501 (2006). 

[35] V. Berezinsky and A. Gazizov, ApJ 643, 8 (2006). 

[36] V. Berezinsky and A. Z. Gazizov (2007), arXiviastro- 
ph/0702102. 

[37] V. Springel and L. Hernquist, MNRAS 339, 289 (2003). 
[38] KASKADE collaboration: K.-H. Kampert et al., Acta 

Physica Polonica B 35, 1799 (2004). 
[39] Akeno collaboration: M. Nagano et al., J. Phys. G 18, 

423 (1992). 

[40] AG AS A collaboration: M. Takeda et al., Phys. Rev. L 

81, 1163 (1998). 
[41] HiRes collaboration: R. U. Abbasi et al., Phys. Rev. L 

92, 151101 (2004). 
[42] Fly's Eyes collaboration: D. J. Bird et al., ApJ 424, 491 

(1994). 

[43] E. Waxman and J. Miralda-Escude, MNRAS 472, 89 
(1996). 

[44] G. Sigl (2007), arXiv:astro-ph/0703403. 



